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Abstract — This study addresses the retrieval of spatially con- 
tiguous canopy cover and height estimates in southwestern US 
forests via inversion of a geometric-optical (GO) model against 
surface bidirectional reflectance factor (BRF) estimates from the 
Multi-angle Imaging SpectroRadiometer (MISR). Model inver- 
sion can provide such maps if good estimates of the background 
bidirectional reflectance distribution function (BRDF) are avail- 
able. The study area is in the Sierra National Forest in the Sierra 
Nevada of California. Tree number density, mean crown radius, 
and fractional cover reference estimates were obtained via analysis 
of QuickBird 0.6 m spatial resolution panchromatic imagery using 
the CANopy Analysis with Panchromatic Imagery (CANAPI) 
algorithm, while RH50, RH75 and RH100 (50%, 75%, and 100% 
energy return) height data were obtained from the NASA Laser 
Vegetation Imaging Sensor (LVIS), a full waveform light detection 
and ranging (lidar) instrument. These canopy parameters were 
used to drive a modified version of the simple GO model (SGM), 
accurately reproducing patterns of MISR 672 nm band surface re- 
flectance (mean RMSE = 0.011, mean R 2 = 0.82, N = 1048). 
Cover and height maps were obtained through model inversion 
against MISR 672 nm reflectance estimates on a 250 m grid. 
The free parameters were tree number density and mean crown 
radius. RMSE values with respect to reference data for the cover 
and height retrievals were 0.05 and 6.65 m, respectively, with R 2 
of 0.54 and 0.49. MISR can thus provide maps of forest cover and 
height in areas of topographic variation although refinements are 
required to improve retrieval precision. 
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I. INTRODUCTION 

T HIS study uses high quality reference data derived from 
high spatial resolution (0.6 m) panchromatic imagery 
and full waveform lidar to investigate the accuracy with which 
forest canopy cover and height can be retrieved by adjust- 
ment of a geometric-optical canopy reflectance model against 
multiangle 672 nm reflectance data from MISR, a moderate 
resolution imager on NASA’s Terra satellite. Forest canopy 
cover and height are important first-order canopy structure pa- 
rameters in basic and applied ecological research. They reflect 
tree density, successional stage, canopy response to stand-re- 
placing fires and other disturbance, forest health, susceptibility 
to fire, species distribution, wildlife habitat, and carbon storage 
in aboveground woody biomass. In the southwestern United 
States the need for forest canopy cover, height, and biomass 
mapping capabilities has never been more pressing, with open 
canopy forests expected to come under increasing stress from 
climate-related disturbance, e.g., unusually severe and exten- 
sive wildfires, mortality from insect and pathogen outbreaks, as 
well human population and management factors (e.g., grazing, 
fire suppression) in the coming decades [l]-[6]. Since the 
mid-1970s the southwestern US has warmed rapidly and pre- 
cipitation has decreased leading to more intense drought [7]; 
moreover, there is a consensus that trends towards warming 
and increased potential evapotranspiration will cause negative 
water balances by the middle of the century [8], [9]. This 
will affect forests importantly: higher elevations — where most 
southwestern forests exist — have experienced greater warming 
than low elevations. The implication is that if temperature and 
aridity follow their expected trajectories, these forests will 
experience substantially reduced growth during this century 
[6]. Warmer winters and earlier melting of the snowpack and 
drying of soils and fuels in the spring are already linked to 
increasingly extensive and severe wildfires in western U.S. 
forests [1]. Some fires are stand-replacing whereas others are 
less damaging: canopy cover and height maps can help in 
remotely determining the severity of a bum. In addition, cover 
and height are important for assessing habitat. For example, 
current management on all Sierra Nevada National Forests 
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uses canopy cover and tree size guidelines based on the need 
to provide habitat for sensitive species such as the California 
spotted owl (Strix occidentalis occidentalis). At present the 
Forest Service often estimates habitat suitability from sparse 
plot data using a model, the Forest Vegetation Simulator [10]. 

Lidar, radar, and passive multiangle imaging all have a role 
to play in mapping these forest parameters at regional scales. 
For example, the Geoscience Laser Altimeter System — a 
large footprint profiling waveform lidar that was flown on the 
Ice, Cloud, and land Elevation Satellite — has provided good 
estimates of canopy height over large areas [11] but did not 
provide contiguous mapping. Airborne instruments such as 
the NASA Laser Vegetation Imaging Sensor (LVIS, [12]) can 
provide robust, precise, and contiguous estimates of canopy 
height and other structural measures [13] but only sporadi- 
cally and over relatively small areas. Radar data can provide 
estimates of canopy height when additional elevation data 
are available [14]; and aboveground biomass when calibrated 
with lidar or ground measurements [15]. Multiangle passive 
solar wavelength reflectance data from moderate resolution 
instruments provides a third option that has some advantages 
but also presents important challenges, discussed below. 

For assessing changes in forest at regional scales and larger 
it is necessary to adopt an approach that will provide informa- 
tion on a regular basis and with as long a record as possible. The 
estimates must have both the accuracy and precision sufficient 
to differentiate real from spurious trends and to resolve changes 
from forest degradation or re-growth as well as large changes 
from harvesting or fire. It is relatively straightforward to de- 
tect total forest loss from a large, severe wildfire with a tech- 
nology that provides a cover estimation precision of ±30% but it 
would be impossible to say anything meaningful about gradual 
or subtle changes, such as those arising from selective logging 
and recovery, especially over periods of less than 10 years. This 
implies that there is utility in data from moderate resolution in- 
struments (i.e., those with instantaneous fields-of-view (IFOVs) 
of 100-300 m) that provide a long record and have good spa- 
tial and temporal coverage. However the data are difficult to 
interpret in terms of structural attributes such as canopy cover 
and height; canopy physical quantities must be inferred from 
variation in the remotely-sensed signal with respect to wave- 
length, polarization, and/or viewing and illumination angles, 
rather than from information available in the spatial domain, as 
with photo grammetric methods. Moreover, the results are diffi- 
cult to validate because of the large ground resolution element of 
>=250 m: it is difficult and expensive to obtain ground-based 
measurements that cover the range of conditions sufficiently 
well. Plot data that may be used with high- and medium-res- 
olution imagery — such as the U.S. Forest Service’s Forest In- 
ventory Analysis (FIA) survey data — are at the wrong scale for 
use with moderate resolution imagery because the 36 m diam- 
eter FIA plot size is almost an order of magnitude too small. 

There are several approaches to mapping canopy cover and 
height over large areas using moderate resolution multiangle 
satellite remote sensing but they may be divided into two broad 
categories: empirical and physical. The empirical approach re- 
lies on calibration of models whose parameters while poten- 
tially meaningful are not physical quantities (e.g., the coeffi- 
cients of a regression equation or the connection weights used 


in an artificial neural network). The physical approach exploits 
models that describe the variation in the remotely-sensed signal 
with respect to wavelength, polarization, and/or viewing and il- 
lumination angles. Both categories involve models of varying 
complexity, from simple to highly complex and with few or 
many parameters. Heiskanen [16] and Kimes et al [17] provide 
examples of empirical approaches to interpretation of multi- 
angle data (exploiting regression and artificial neural networks, 
respectively); while Widlowski et al [18], Schull et al [19], 
Zeng & Schaepman [20], Chopping et al [21]-[23], Wang et al 
[24], and Laurent et al [25] are some examples of physical ap- 
proaches (BRDF models, geometric-optical models, directional 
escape probabilities). These studies have shown that it is pos- 
sible to access the information on canopy structure that is en- 
capsulated in multiangle reflectance data. 

The approach pursued in this study involves adjusting a geo- 
metric-optical (GO) model in a framework in which the con- 
tribution of the (non-tree) background at different viewing and 
illumination angles is estimated for each location using a dy- 
namically-estimated bidirectional reflectance distribution func- 
tion (BRDF), rather than a static surface BRDF or brightness 
value [2 1 ]— [23]. A BRDF model provides both reflectance mag- 
nitude and anisotropy (i.e., the angular distribution of the out- 
going radiation field with respect to the surface and the solar 
direction), rather than only the isotropic, or diffuse, scattering 
component. Specification of the background BRDF — or at least 
the part of it corresponding to the viewing and illumination an- 
gles of the observations — is critical for accurate GO canopy re- 
flectance modeling. The background has long been recognized 
as an important component of the remote sensing problem [26], 
[27] and remains an area of active research [28]. In this study 
the background BRDF is estimated via least-squares regression, 
using the kernel weights of a BRDF model that has been ad- 
justed against the available set of bidirectional reflectance factor 
(BRF) observations. 

Two previous published studies have demonstrated the 
utility of red band multiangle BRF data from MISR in such a 
GO model inversion framework (hereafter termed MISR/GO) 
for forest canopy mapping. In a first study over an area of 
~200,000 km 2 in New Mexico and Arizona, retrievals of 
fractional cover (fcov ) and mean canopy height were assessed 
against United States Department of Agriculture, Forest Service 
Interior West (FS-IW) map series [21]. However the FS-IW 
estimates were derived using an empirical modeling technique 
that ingested moderate resolution multispectral products from 
the NASA MODerate resolution Imaging Spectroradiometer 
(MODIS), as well as other geospatial inputs and so are not 
suitable for validation purposes, especially for structural pa- 
rameters such as canopy height [22], [29]. A second study 
assessed canopy height retrievals against high resolution dis- 
crete return lidar heights from the NASA Cold Land Processes 
Experiment (CLPX) in the Colorado Rockies. Forest height 
RMSE distributions with respect to the CLPX lidar estimates 
were centered between 2.5 and 3.7 m while R 2 distributions 
were centered between 0.4 and 0.7 (N = 57) [22]. 

While good compatibility was found between MISR/GO 
cover, height, and aboveground biomass retrievals and the 
FS-IW estimates and height retrievals were accurate with 
respect to high resolution lidar height measurements, initial 
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TABLE I 
Data Sets Used 


Date 1 

Extent (ha) 

Description 

- 

652,406 2 * 

MISR Ancillary Geographic Product 

08/11/08 

652,406 2 

MISR 1.1 km LAND Product 

08/11/08 

652,406 2 

MISR Geometric Parameters Product 

08/11/08 

652,406 23 

MISR Terrain- projected Spectral Radiance 

09/26/08 

6,536 4 

LVIS Canopy and Ground Elevation 

06/25/03 

7,380 5 

QuickBird 0.6m Panchromatic Image 


Filename(s) 

MISR_AM 1_AGP_P04 1_F0 l_24.b60-61.hdf 

MISR_AMl_AS_LAND_P041_0046009_F07_0022.b60-61.hdf 

MISR_AMl_GP_GMP_P041_0046009_F03_0013.b60-61.hdf 

MISR_AM 1_GRP_TERRAIN_GM_P04 l_0046009_nn_F03_0024 .b60-6 1 .hdf 

LVIS_US_CA_day4_2008_VECT_2008 1 120.lge 

031UN25183856-P2AS_R2C1-052222275010_01_P001.TIF 


1 Dale of acquisition by the remote sensing instrument 

: MISR blocks 60-61 were mapped into the area bounded by 37° 10' 34"N. 1 19° 17 4<TW by 36° 31’ 26"N. 1 18° 18' 22"W. 

5 nn indicates the codes for the nine MISR cameras (AN. A A. AF. BA. BF. CA. CF. DA DF). 

4 Irregular area sampled; the validation area corresponded to lire area of the QuickBird subset minus the NE and SW comers wliere no LVIS data were acquired. 

5 The area is bounded by 37° 01' 18"N. 1 19° 09* 09 'W by 36° 57 00 ,r N. 1 19° 02' 57"W. 


attempts to map the entire southwestern US using a MISR/GO 
method produced anomalous results in some places. For ex- 
ample, in locations where the topography is complex and 
sometimes extreme, cover was severely overestimated and 
height severely underestimated [30]. This can happen because 
of the non-uniqueness of solutions when adjusting a model 
using numerical methods and the pattern of modeled BRFs 
is very similar with different parameter sets (i.e., the same 
minimum error within the defined convergence precision can 
be obtained with more than one solution). Such artifacts may be 
owing to model inadequacy, poor estimates of the background 
BRDF, geometric effects of topography, residual atmospheric 
or intrinsic noise in the BRF data, and/or problems with the 
inversion protocol (e.g., sensitivity to the starting point). 

This study uses robust high quality reference data derived 
from analysis of high spatial resolution (0.6 m) panchromatic 
imagery and full waveform lidar to investigate the precision and 
accuracy with which forest cover and height can be retrieved by 
adjustment of a GO model against multiangle reflectance data 
from MISR. Lidar data acquired from the air provide the best 
means of adequately assessing MISR/GO height retrievals; in 
this study height metrics from the NASA LVIS full waveform 
lidar with a nominal spot size of 20 m were used. Accurate es- 
timates of tree density and crown cover, as well as first order 
tree heights from shadowing were provided by the CANopy 
Analysis with Panchromatic Imagery (CANAPI) algorithm that 
was developed to provide maps of these parameters using high 
spatial resolution (0.5-1 m) panchromatic imagery [31]. These 
data are required for validation of MISR/GO retrievals and also 
to drive forward-modeling runs that provide an indication of 
whether the model is able to accurately reproduce BRF patterns 
at MISR observation geometries; they are described more fully 
below. 

II. Methods 

In this study tree number density, mean crown radius de- 
rived from QuickBird imagery with CANAPI and canopy height 
from analysis of LVIS waveforms were used to drive a geo- 
metric-optical canopy reflectance model in forward mode (to 
test the ability of a model to predict the 672 nm wavelength ra- 
diation field above the canopy as seen by MISR) — as well as 



Fig. 1 . The location of the study area in the Sierra Nevada National Forest in 
California, USA (a) over a crown cover map (b) MISR multi-angle composite 
(RGB = vol, geo, iso kernel weights). The irregular area indicates the area sur- 
veyed by LVIS on day 4 of the September 2008 campaign. The small rectangle 
shows the extent of the area used to collect high quality reference data. Black 
indicates data missing owing to cloud. 

in inverse mode (to extract background BRDFs for calibration 
of the regression models that provide the weights required for 
calculation of the background contribution for each location, 
prior to adjusting the model to obtain fractional crown cover 
and mean canopy height). The following sections describe the 
data and procedures followed. Please see Table I for a list of the 
data sets used and their spatio-temporal extents. 

A. Study Area 

The study area is centered on the area surveyed by LVIS in 
September 2008 in the Sierra National Forest in the southern 
Sierra Nevada Mountains of California, USA (Fig. 1). It encom- 
passes an area of 65,625 ha (37° 10' 34" N, 119° 17' 40" W by 
36° 31' 26" N, 118° 18' 22" W) with elevation ranging from 
~850 to 2,700 m. The area includes upland forest in the western 
part and high desert grassland and scrub in the eastern part. 
The upland area covers a wide range of vegetation associations 
including fir, pine, mixed conifer, mixed hardwood and conifer 
forests sparsely interspersed with meadows. Common tree 
species include red fir ( Abies magnifica ), white fir ( Abies con- 
color ■), ponderosa pine ( Pinus ponderosa ), Jeffrey pine ( Pinus 
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jeffreyi ), and incense cedar ( Calocedrus decurrens ) [32]. Tree 
canopy cover can range from completely open in meadows, 
rocky ridges, and scree slopes to very dense (but rarely ex- 
ceeding 50%, or fractional cover = 0.5). The topography is 
highly variable, often complex, and sometimes extreme, with 
many steep slopes, some well over 40° . There is high variation 
in slope and aspect within 250x250 m mapped MISR cells, 
with standard deviations between 0°-13° (slope) and 2°-158° 
(aspect). On the mountain range precipitation ranges from 
500-2030 mm during fall, winter, and spring and occurs mostly 
as snow above 1828 m. Summers are dry with low humidity, 
and temperature averages 5.5-1 5.5° C with a growing season 
lasting 20-230 days. On the foothills precipitation ranges 
from 510-1020 mm, temperature averages 13-18°C, and the 
growing season lasts 200 to 320 days [33]. 


B. MISR Data Sets and Processing 

The MISR instrument is flown on Terra, the first of NASA’s 
Earth Observing System satellites that was launched in De- 
cember 1999. The instrument consists of nine pushbroom 
cameras arranged to view along-track that acquire image data 
with nominal view zenith angles relative to the surface refer- 
ence ellipsoid of 0.0°, ±26.1°, ±45.6°, ±60.0°, and ±70.5° 
(forward and aft of the Terra satellite) in four spectral bands 
(446, 558, 672, and 866 nm) at 1.1 km spatial resolution. The 
672 nm (red) band images are also acquired with a nominal 
maximum cross-track ground spatial resolution of 275 m in all 
nine cameras and all bands are acquired at this resolution in the 
nadir camera [34]. Except as otherwise noted, MISR data prod- 
ucts are obtained using the MISR Order and Customization Tool 
hosted at the NASA Langley Atmospheric Science Data Center 
(http://10dup05.larc.nasa.gov/MISR/cgi-bin/MISR/main.cgi). 

The high resolution (275 m) surface BRF is obtained by first 
performing a linear regression of the red band 1.1 km MISR 
surface BRFs against the red band 1.1 km MISR top-of-atmos- 
phere BRFs to obtain the regression coefficients for each of the 
nine camera view data sets. The regression is performed on the 
data within each 17.6 km region, producing a grid of regression 
coefficients, which then are smoothed to eliminate discontinu- 
ities between regions. The smoothed regression coefficients are 
then applied to the red band 275 m top-of-atmosphere MISR 
BRFs to produce red band 275 m land surface BRFs at the nine 
camera angles. Top-of-atmosphere BRFs contaminated by sun 
glint from water surfaces are excluded from input to the regres- 
sion coefficient calculations. 

The surface BRF regression algorithm utilizes four MISR 
data products: the MISR Level- 1B2 Terrain-projected Radiance 
product (MI1B2T); the MISR Level-1B2 Geometric Parame- 
ters product (MI1B2GEOP); the MISR Level-2 Land Surface 
product (MIL2ASLS); and the Ancillary Geographic product 
(MIANCAGP). The MI1B2T product provides 275 m terrain- 
projected top-of-atmosphere radiance in all 4 spectral bands for 
the nadir view; and in the red band for off-nadir views. Sun and 
view geometries from the MI1B2GEOP product are combined 
with a land/water mask from the MIANCAGP product to de- 
termine view angles and surfaces susceptible to sun glint. The 
MIL2ASLS product provides surface BRFs at 1.1 km. 
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Fig. 2. Schematic showing the principles on which the CANAPI algorithm iso- 
lates tree or shrub crowns, finds crown radii, and estimates heights from shadow 
lengths. The algorithm is implemented in the Image J image processing package 
from the U.S. National Institutes of Health. 


The MISR data used in this study were from Terra orbit 
46009 (August 1 1 , 2008), path 4 1 , blocks 60-6 1 . The data were 
resampled onto a 250 m grid in the Universal Transverse Mer- 
cator projection, Zone 11N, WGS84 spheroid/datum. The latest 
operational MIL2ASLS product (version 22) available from the 
Langley data center lacked adequate coverage in the study area, 
primarily due to a topographic complexity constraint imposed 
by the MISR aerosol algorithm ([35], sections 3. 3. 1.2.2 and 
3. 3. 8. 2. 5). Relaxing the regional and subregional topographic 
complexity thresholds to 1000 m (from 500 m and 250 m, 
respectively), substantially improved coverage in mountainous 
regions. Additional coverage was gained by allowing the 
aerosol retrieval to proceed without the ±70.5° views, when 
doing so improved coverage. This latter modification is also 
relevant in mountainous regions due to topographic obscura- 
tion of the oblique camera views. The improved MIL2ASLS 
product was used as input to the surface BRF regression. 

C. Forest Canopy Cover and Statistics From Panchromatic 
Imagery 

The CANopy Analysis with Panchromatic Imagery 
(CANAPI) algorithm allows the calculation of estimates 
of crown cover, mean crown radius, tree number density, 
background (non-tree) brightness, and tree height, over large 
areas [31]. The algorithm identifies crowns in high resolution 
panchromatic imagery by finding the crescent-shaped areas 
that correspond to illuminated crowns. It accomplishes this by 
rotating the original image so that the sun direction is at the 
top and applying normalization and convolution operations 
to emphasize sunlit crown pixels (Fig. 2). CANAPI produces 
maps of idealized (circular) crowns, crown centers, radii, and 
tree height from shadow-following, wherever possible (Fig. 3). 
Statistics such as tree number density and mean crown radius 
can then be calculated for a window of user-selected dimen- 
sions, with the mean radius calculated after removing a few 
values at the edge of the distribution because the maximum is 
often the result of misidentification of several crowns as one 
and the minimum is constrained by pixel size. 
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Fig. 3. Typical CANAPI crown (circle) and tree shadow (line) detections over 
QuickBird 0.6 m panchromatic images in the Teakettle Experimental Forest, 
Sierra National Forest, California. Shadows that are truncated by tree crowns or 
the edge of the image are not used in tree height calculation. The imagery was 
acquired June 25, 2003. 


Here, CANAPI estimates of number density and mean crown 
radius were calculated for 250x250 m windows corresponding 
to mapped MISR cells, using QuickBird 0.6 m spatial resolution 
panchromatic imagery acquired on June 25, 2003. The area cor- 
responds to 32x34 MISR cells, providing a reference data set 
with N = 1088. This particular image was selected because it 
is one of the few available that provides extensive views of the 
surface during the summer season. There is the possibility that 
the forest may have changed over the five-year period between 
the acquisition of the QuickBird imagery and the other data sets, 
owing to disturbance. In Sierra Nevada mixed conifer forest 
the main disturbances affecting canopy conditions are fire and 
bark beetle mortality; wind damage, pests and snow breakage 
are very rare, localized events. While it is not absolutely cer- 
tain that there were no major disturbance events in the valida- 
tion area over the 2003-2008 period, this seems unlikely: it is 
known that there were no fires in the Teakettle Experimental 
Forest, a 1300 ha old growth watershed about 2 km to the east 
that has been intensively studied (e.g., [36]). Bark beetles and 
other pests and pathogens were studied at Teakettle over this pe- 
riod [37]: while there were chronic levels of bark beetle damage 
these affected a very small fraction of trees and there were no 
significant droughts — which are a catalyst for increased beetle 
mortality — during this period. In the absence of these major dis- 
turbances there were likely only small changes in primary forest 
canopy structural properties (crown cover, mean height) over 
the period June 2003-September 2008. 

CANAPI provides very accurate estimates of fractional 
cover. When these estimates were compared with those derived 
from the extensive and detailed Teakettle Ecosystem Experi- 
ment (TEE) database of field measurements collected between 


1998 and 2001 for fifteen 200 m 2 plots [36], they yielded an 
absolute RMSE of 0.03 and an adjusted R 2 of 0.92, significant 
at the 99% level [31]. Maximum CAN API-derived heights 
showed a good spatial agreement with LVIS RH100 (100% 
energy return) maximum canopy height estimates interpolated 
onto a 20 m grid, although the CANAPI map includes far fewer 
very tall trees and the relationship was rather weak over the 
entire range of RH100 values, with obvious underestimation 
and anomalies at low and high values owing to the inability to 
obtain heights for all trees falling in the calculation window 
(shadows may be truncated by crowns or the edge of the 
window). However when the range 3 m to 60 m (61% of the 
total) was considered then the R 2 was 0.94, significant at the 
99% level, with a relative RMSE of 0.25 and an absolute RMSE 
of 13.9 m that reflects the divergence from the 1:1 line [31]. 

D. Lidar Canopy Heights 

The lidar data used in this study were collected by the Laser 
Vegetation Imaging Sensor (LVIS) [12]. LVIS is a full-wave- 
form laser altimeter system optimized to measure canopy 
structure across large areas. In September of 2008, as a part of 
an experiment for NASA’s Deformation, Ecosystem Structure 
and Dynamics of Ice (DESDynl) mission, LVIS mapped a 
large portion of the Sierra National Forest with a 1-2 km wide 
swath in a series of parallel flight tracks (Fig. 1). Flying at high 
altitude LVIS produced millions of (nominally 20-m diameter) 
footprints at the surface. The data from the LVIS surveys used 
here were collected on Day 4 of the campaign (September 
26, 2008) and cover approximately 7,400 ha (37° 01' 16"N, 
119° 09' 09" W by 36° 56' 60"N, 119° 02' 58" W). They are 
entirely within the area covered by the MISR imagery. The 
waveforms were analyzed by the LVIS science team to provide 
ground elevation and canopy height estimates and made avail- 
able via the LVIS web server at NASA Goddard Space Flight 
Center [38]. The products used were the ground elevation 
estimates (“LGE”) with respect to the International Terrestrial 
Reference Frame 2000 (horizontal coordinates)/WGS84 El- 
lipsoid (vertical coordinates) and the canopy relative height 
metrics RH50, RH75, and RH100 (“LCE”). In the LVIS 
product, the height at which the indicated percentage of energy 
has been returned is estimated; thus RH100 is a measure of 
maximum tree height within the lidar footprint while RH50 is 
the height of median energy (HOME). RH75 is the height at 
which 75% of the energy is returned to the sensor. The mean 
incidence angle of the LVIS data set was 2.8° with a standard 
deviation 1 .4° ; none of the data were acquired at an incidence 
angle greater than 10°. The LVIS RH50, RH75, and RH100 
values within each mapped MISR 250x250 m cell were av- 
eraged to obtain first-order canopy height reference data sets, 
thus these data are not original LVIS data but quantities derived 
from them. The LGE elevations were used to construct a 37 m 
horizontal resolution digital elevation model (DEM, Fig. 4) that 
was used to effect first-order corrections to the RH100 values 
for slope. 

E. Modified Simple Geometric-Optical Model 

Geometric-optical (GO) models are able to resolve statistical 
distributions of discrete objects within an instrument’s instan- 
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Fig. 4. Topography of the Sierra National Forest study site (a) digital elevation 
model (DEM) using linear interpolation on LVIS ground returns (b) slope map 
generated from the DEM (c) aspect map generated from the DEM. Only valid 
data in the center of the maps — as shown by the solid lines — were used. 


taneous field-of-view (IFOV) [39]. Simple geometric-optical 
models treat the surface as an assemblage of discrete objects 
of equal radius, shape and height, evenly distributed within a 
spatial unit. A tree or shrub crown is represented by a spheroid 
whose center is located at a specified mean height above a 
background with defined reflectance magnitude and anisotropy 
(i.e., the background is represented by a BRDF model). These 
models predict the top-of-canopy reflectance response to im- 
portant canopy physical parameters (e.g., plant number density, 
foliage volume, mean canopy crown height, radius, and crown 
shape, background brightness and anisotropy) as a linear 
combination of the contributions from sunlit and viewed, and 
shaded and viewed crown and background components [40], 
[41]: 

R = G.k G + C.k c + T.k T + Z.k z (1) 

where R is bidirectional spectral reflectance factor; k G ,k G , 
and k z are the GO modeled proportions of sunlit background, 
sunlit crown, shaded crown and shaded background, respec- 
tively; and G , C , T, and Z are the contributions of the sunlit 
background, sunlit crown, shaded crown, and shaded back- 
ground, respectively. GO functions have been used in Li-Ross 
BRDF models [42], [43] and are particularly appropriate for 
the exploitation of solar wavelength remote sensing data ac- 
quired at differing viewing and/or illumination angles (i.e., 
multiangle remote sensing data) because the proportions of 
sunlit and shaded crown and background in the remote sensing 
instrument ground-projected IFOV vary with both viewing and 
illumination angles and canopy configuration. 

GO models calculate the proportions of sunlit background 
(k G ) and crown (k G ) with respect to illumination and viewing 
directions. These are calculated exactly for the principal and 


perpendicular planes and approximated away from these, as in 
(2) and (3), respectively: 

g — A7rr 2 {sec ^+sec $' v — 0(i^i ,<£>)} ^2^ 

k c = (l _e- Wsec <) 1(1 + cose') (3) 

where # and ip indicate zenith and azimuth angles where the 
subscripts i and v indicate illumination and viewing angles; A is 
the tree number density; r is the average radius of tree crowns; 
e' is the transformed scattering phase angle: 

cos s' = cos cos + sin sin $' v cos tp (4) 

and O is the overlap area between the shadows of illumination 
and viewing [43]: 

O = 1/7t(£ — sin t cos t) (sec + sec $' v ) (5) 

where t is a parameter that indirectly expresses the locations of 
the end points of the line that intersects the shadows of viewing 
and illumination. This allows k G to be expressed in a way that 
depends only on the value of t (see [43]). The prime indicates 
equivalent zenith angles obtained by a vertical scale transfor- 
mation based on b/r, in order to treat spheroids as spheres [43]: 

= tan -1 (6/r tan$) (6) 

These functions thus depend on the parameters b/r (vertical 
crown radius/horizontal crown radius) and hjb (height of 
crown center/vertical crown radius) which describe the shape 
and height of the crown. The model’s primary structural param- 
eters are therefore: tree number density (A), mean crown radius 
(r), crown vertical to horizontal radius ratio (b/r), and crown 
center height to vertical radius ratio (h/b). The crown signature 
may be represented by a volume scattering function [44], in 
which case an in-crown LAI or foliage volume parameter may 
also be specified. 

In this study a modified version of the simple geometric 
model (SGM), a GO model incorporating a dynamic back- 
ground (calculated for each location separately), was used 
[21]— [23]. It is formulated as 

R = GRTLs{$i, $v, ip).k G (iH, tiv, if) 

+C crown .k c {'&i, tiv, ip) + S TZ .k TZ ('di , Av, f) (7) 

where $i, and ip are the view zenith, solar zenith and 
relative azimuth angles, respectively; k G , k G and k Z r are the 
calculated proportions of sunlit and viewed background and 
crown, and shaded components (shaded ground and crown), 
respectively; Grtls is the background contribution repre- 
sented by the RossThick-LiSparse-reciprocal (RTLS) BRDF 
model; C crown is the red wavelength bidirectional reflectance 
signature for crowns (set at 0.03 in this study); and Srz is 
the estimated contribution of shaded components (here set to 
0.02). Note that the crown signature differs from the spectral 
reflectance of fir leaves (~0.05; [45]) because the observed 
crown area in the instrument IFOV is not composed simply 
of leaves but includes non-negligible in-crown shadows. The 
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crown and shadow signature values are estimated but they are 
not arbitrary. The model’s physical upper canopy parameters 
are tree number density (A, trees per square meter), mean 
crown radius (r), crown vertical to horizontal radius ratio 
(b/r), crown center height to vertical radius ratio (h/b), and 
crown spectral reflectance. The model used in this study differs 
from the previous version in several ways: 

• The shaded components T and Z (( 1 )) that were previously 
assumed black (as in the kernel-driven BRDF models) are 
no longer discarded because in these tall forests shaded 
ground and crown can account for a relatively large pro- 
portion of a nadir- viewing instrument’s IFOV (Fig. 3(c)); 
instead they are allocated a signature, S T z , of 0.02; 

• The 3 -parameter RTLS BRDF model replaced the 4-pa- 
rameter modified Walthall model [46], [47] in order to 
make assessment of extracted and predicted backgrounds 
more straightforward; and to increase the robustness of 
background BRDF retrievals (with four kernels it is more 
likely that their contributions will be confounded); 

• The Ross volume scattering function with its in-crown LAI 
parameter is replaced by a crown red wavelength signature 
of 0.03, based on a leaf red wavelength reflectance value 
of 0.05 for fir [45] minus a small decrement to account for 
in-crown shadows. 

F. Isolation of the Background Contribution 

The “background” includes everything within the mapped 
IFOV of the instrument that is not part of the upper canopy and is 
here typically composed of many elements, possibly including 
bare soil, rock faces, scree, gravel, crusts, mosses, lichens, and 
understory plants such as grasses, forbs, and sub-shrubs (Fig. 5). 
In order to estimate background brightness and anisotropy in 
the various MISR illumination/viewing configurations, linear 
multiple regression based on a number of calibration sites was 
used. The independent variables are RossthiN-LiSparse-recip- 
rocal (RNLS) BRDF model red band isotropic (iso), geometric 
(geo), and volume scattering (vol) kernel weights, obtained by 
adjusting the model against MISR red band BRFs in all nine 
cameras using the Algorithm for Modeling Bidirectional Re- 
flectance Anisotropies of the Land Surface (AMBRALS) code 
[48], with the objective function Min(|RMSE|). Since it is dif- 
ficult to find 250x250 m areas that are partly vegetated but do 
not contain shrubs or trees, background contributions must usu- 
ally be extracted: provided with a set of MISR red band BRFs, 
estimates of upper canopy statistics from CANAPI (and op- 
tionally, LVIS), and the GO model, an optimization algorithm 
can extract the best-matching background for each site. In other 
words, when the contribution of the upper canopy in the ob- 
served MISR BRF pattern has been accounted for, the remaining 
contribution must be from the background (Fig. 6). Background 
BRDFs retrieved in this way are hereafter referred to as “op- 
timal” backgrounds. 

Two methods were used to find sets of background RTLS 
model red band kernel weights (i.e., the dependent variables 
in the regression equations). First, an automated approach was 
followed: since CANAPI provides canopy statistics for large 
areas it is possible to extract kernel weight sets for many con- 
tiguous 250 x 250 m windows corresponding to MISR grid cells. 



Fig. 5. A selection of 250 x 250 m image chips from Google Earth showing 
the wide range of non-tree- backgrounds for which the GO modeling frame- 
work must attempt to predict BRDF: (a) rocky/stony scrub; (b) meadow with 
lichens; (c) senescent sparse grasses with rocks; (d) lush meadow; (e) senescent 
grasses; (f) senescent grasses with exposed soil; (g) low shrubs, rocks/stones, 
exposed soil, grasses; (h) highly stratified uplifted and tilted rock face; (i) senes- 
cent grasses with exposed bright mineral soil; (j) bright exposed caliche with 
rock outcrops and stones; (k) wet meadow grasses and low shrubs; (1) talus 
slope. 

Here 1048 windows from the part of a QuickBird scene from 
June 25, 2003 that included LVIS canopy height metrics were 
used (the “validation area”). To find the optimal backgrounds, 
for each location the GO model A and r parameters were set 
to the CANAPI-derived values, h/b was fixed at 1.0 (typical 
for conifers), and b/r was adjusted so that mean crown height 
(h + b) matched the corresponding LVIS RH50 value. RH50 
was used because it is more closely related to GO model mean 
canopy height than RH75 or RH 100 [13], [49]. The RTLS back- 
ground BRDF model kernel weights were then adjusted to min- 
imize the absolute RMSE between the GO model and MISR 
BRFs. In the second approach, an extensive visual scan of high 
resolution imagery was performed in Google Earth for the area 
of the MISR scene, seeking 250x250 m cells with very low 
(generally <10%) canopy cover. The same extraction proce- 
dure was followed, though mean crown radius was estimated 
from the imagery, tree number density set to give a match with 
fractional cover (measured via thresholding), and b/r was set to 
1.0 and 4.0 for broadleaf and coniferous trees, respectively. 

Having obtained the best-fitting RTLS background model 
kernel weights with respect to the MISR data, regression 
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coefficients were obtained for each of the model parameters 
using the RNLS kernel weights as the independent variables. 
The resulting regression equations can be used to obtain esti- 
mates of the background response prior to adjustment of the 
GO model [21]— [23]. For the automated approach, the set of 
locations used to obtain the regression coefficients was reduced 
by considering only those locations where the CANAPI crown 
cover estimate was <15%. 

G. Model Evaluation in Forward Mode 

In forward mode the GO model was driven to predict red band 
BRFs for all available MISR cameras (in these data ranging 
from six to nine looks but usually eight or nine) for all loca- 
tions with reference data, using tree number density and mean 
crown radius from CANAPI; mean b/r from adjustment against 
LVIS RH100 (because modeled height depends on b/r); and the 
RTLS BRDF model background (predicted by red band RNLS 
BRDF model kernel weights). This was accomplished for 1048 
of the 1088 250x250 m areas within the QuickBird sub-scene 
as complete LVIS data were not available for 40 MISR pixels. 
The resulting BRFs were then compared with the observed pat- 
terns from MISR. 

H. Model Inversion Protocol 

Retrieval of fractional cover and mean canopy height for 
each mapped grid location was performed by adjusting the 
model against the MISR data using a numerical optimization 
algorithm. Forward differencing was used for estimates of par- 
tial derivatives of the objective function and a Newton search 
method was used at each iteration to decide which direction 
to pursue in the parameter space [50]. The objective function 



Fig. 7. Flowchart of the operations performed to retrieve and validate fractional 
crown cover and mean canopy height via GO model inversion against MISR 
672 nm band BRFs. 


was Min (| RMSE | ) with no weighting of the error terms with 
respect to viewing angle. No constraints were imposed. Tree 
number density (A) and mean radius (r) were left as adjustable 
parameters, with hjb and b/r set to 1.0 (typical for conifers) 
and 4.0 (tall trees), respectively. Retrievals were thus effec- 
tively for fractional crown cover (a function of A and r) and 
mean canopy height (a function of r, since model height is 
controlled by the terms hjb and b/r that are both fixed). MISR 
red band data in all available cameras were used. The number 
available is usually nine but is occasionally six, seven, or eight. 

The initial values of A and r were set to 0.004 (250 trees per 
MISR grid cell of 62,500 m 2 ) and 3.0 m, respectively. Note that 
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Fig. 8. (a)-(c) Predicted versus optimal background BRDF model red band isotropic, volume, and geometric scattering kernel weights for the 28 background 
calibration sites; (d)-(f) for backgrounds predicted using kernel weights plus the green and near- infrared BRFs from the MISR nadir camera. 


effective values are retrieved: the parameters are internally cou- 
pled and the same BRF pattern is obtained if these two param- 
eters are varied with fractional cover maintained, so it is not 
feasible to adjust both parameters simultaneously to obtain real 
values. However r also determines mean top-of-canopy height 
(h+b): since b/r and hjb are fixed, and A and r are adjusted, the 
inversion effectively provides an estimate of b (vertical crown 
radius) allowing the calculation of an estimate of h, via h/bxb , 
noting that h refers to the mean height of crown centers. Frac- 
tional crown cover ( fcov ) was calculated assuming overlap be- 
tween crowns, according to the Beer-Lambert law (e.g., [10]): 

fcov = 1.0 — exp( — A x 7r r 2 ) (8) 

The retrieved values of A and r were read into raster image 
layers for mapping and analysis, along with the model-fitting 
RMSE, R 2 , fractional cover, and top-of-canopy height. 

The procedures by which the data and GO model are used 
to extract background BRDFs and to invert the GO model and 
validate retrievals are depicted in a flow-chart (Fig. 7). 

III. Results and Discussion 

CANAPI crown cover within the validation area ranges from 
0.0-0. 4 and it is therefore clearly critical that the background 
contribution at the instrument illumination and observation an- 
gles is predicted as accurately as possible. This is a challenging 
task using simple regression methods in view of the diversity 


of background types and each candidate set of coefficients must 
be assessed in both forward and inverse modes. Attempts to ob- 
tain a regression model using the automated approach in which 
all 1048 MISR grid cells for which CANAPI A and r and LVIS 
RH50 data were available did not produce reasonable inversion 
results, even when the CANAPI-derived canopy cover was used 
to restrict the set of locations considered to <15%. There was 
only a small improvement when combining a low cover set ob- 
tained in this way with data from a few individually-selected 
low cover locations outside the validation area. This may be be- 
cause the 1048 cells do not cover a sufficiently wide range of 
background conditions. 

Using 28 grid cells selected for low woody plant cover drawn 
from a larger area provided a much wider variety of background 
types, as indicated in Fig. 5. The predicted background BRDF 
model kernel weights rendered R 2 of 0.99, 0.61, and 0.91 for 
the isotropic, volume-scattering, and geometric kernel weights, 
respectively, with respect to the optimal weights (Fig. 8(a)-(c)). 
This is surprisingly good in view of the wide diversity of 
components that may include stony scrub, lush meadow, 
exposed soil, low shrubs, exposed soil, uplifted rock faces, 
rocks of various sizes, caliche, and talus slopes. Only slightly 
improved results were obtained using kernel weights plus green 
and near- infrared band BRFs from the MISR nadir camera 
(Fig. 8(d)— (f)). The regression results show that not all indepen- 
dent variables contributed to accurate results, for example, the 
RNLS geometric scattering kernel weight was not significant 
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TABLE II 

Background Kernel Weight Coefficients 


Kernel Weight Regression Coefficients 
(bold = not significant @ p=0.05) 


Background model (RTLS-Reciprocal) 

Observed (RNLS-Reciprocal) 

Isotropic 

Volume 

Geometric 

Intercept 

0.008 

0.133 

-0.013 

Isotropic kernel weight 

0.955 

-0.960 

0.112 

Volume scattering kernel weight 

-1.096 

2.879 

-0.489 

Geometric kernel weight 

0.334 

3.485 

0.641 




geometric kernel weight 

Fig. 9. Maps of the predicted background BRDF model kernel weights: 
(a) isotropic (diffuse) scattering; (b) volume scattering; (c) geometric 
scattering. 


in predicting the RTLS background model isotropic kernel 
weight (Table II). Maps of the predicted RTLS background 
BRDF model kernel weights show reasonable distributions, 
with higher isotropic kernel weight values for the snow-cov- 
ered mountain tops in the east-central part of the scene and 
lower values in densely- vegetated areas (Fig. 9). The geometric 
kernel weight also shows a strong spatial correspondence with 
snow cover and is higher over sparse vegetation and lower over 
densely- vegetated areas, while the volume scattering kernel 
weight is relatively high over dense vegetation, as expected, 
although its behavior varies over snow-covered areas. 

Comparison of the model-predicted patterns with those from 
MISR for 1048 locations provided distributions in which 96.6% 


haveR 2 > 0.7 and 97.0% have RMSE < 0.02 [Fig. 10(a), (b)]. 
When the RMSE and R 2 values are plotted against fractional 
crown cover it can be seen that there is only a small change 
in accuracy with cover [Fig. 10(c), (d)], with a mean RMSE 
of 0.011 and a mean R 2 of 0.82. These results suggest that the 
GO model is able to predict the sequence of BRFs observed by 
MISR in all nine cameras rather well, implying that the back- 
ground prediction is effective. Note however that even if the 
background prediction were perfect this is no guarantee of ac- 
curacy on model inversion: there is still a potentially strong de- 
pendence on the optimization algorithm, its scaling and preci- 
sion settings, and the starting point. The distributions of tree 
number density (effective), mean crown radius (m; effective), 
fractional cover (unitless) and mean canopy height (m) are all 
approximately normal (Fig. 11). 

Model adjustment against the MISR red band BRFS resulted 
in mean (standard deviation) model-fitting RMSE and coeffi- 
cient of determination values of 0.009 (0.009) and 0.78 (0.199) 
respectively, with mode values of 0.008 and 0.925 (Table III). 
The inversion algorithm was thus able to fit the model to the data 
well. The MISR/GO fractional cover retrievals showed a mod- 
erately strong relationship to the CANAPI estimates, with an 
R 2 of 0.56 [Fig. 12(a)], significant at the 99% level. The RMSE 
with respect to the CANAPI data is 0.05 over a range that spans 
0.01-0.39. The MISR/GO mean canopy height retrievals show 
a somewhat less strong relationship to the LVIS RH100, RH75, 
and RH50 data, with R 2 of 0.49, 0.38, and 0.44, respectively 
[Fig. 12(b)— (d)], all significant at the 99% level. The RMSE 
with respect to LVIS RH100 is 6.65 meters, considerably higher 
than that obtained under the less strenuous topographic condi- 
tions and lower canopies of the Rocky Mountains in Colorado 
[22], although in this Sierra Nevada landscape the mean and 
range are also much larger: 30.1 m, spanning 3.4-62.0 m. There 
is marked overestimation at low values and underestimation at 
high values with respect to RH100, noting that the RH100 value 
indicates the mean of the set of RH1 00 values within the mapped 
MISR cells. This bias is less marked in the relationship with 
RH75 and is almost negligible in the relationship with RH50 
[Fig. 12(c) and (d), respectively], which may indicate that it 
is not reasonable to compare MISR/GO mean canopy heights 
with the mean of many lidar maximum canopy heights (RH 1 00), 
even though RH100 provided the best match. 
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Fig. 10. Results of the GO model predictions when driven by CANAPI tree number density and mean crown radius, LVIS RH100, and the RTLS BRDF model 
representing the background, for 1048 contiguous locations: (a) histogram of modeled versus MISR RMSE values; (b) histogram of modeled versus MISR coef- 
ficient of determination (R 2 ) values; (c) distribution of RMSE values with respect to CANAPI/QuickBird fractional cover estimates (fcov ); (d) distribution of R 2 
values with respect to CANAPI/QuickBird fractional cover estimates. 



0.06 0.11 0.16 0.21 0.26 0.31 0.36 0.41 19.5 23.3 27.0 30.8 34.5 38.3 42.0 


Fractional Cover Mean Canopy Height (m) 

Fig. 11. Histograms for tree number, mean crown radius, fractional cover, and 
mean canopy height obtained by inversion of the GO model against MISR red 
band BRFs in all nine cameras. 

Representations of typical 250 m grid interval MISR/GO 
maps of model-fitting RMSE, fractional cover, and mean 
canopy height are given in Fig. 13, together with a standard 


TABLE III 

SGM Inversion Statistics 



Model-Fitting 

RMSE 

Model-Fitting 

R 2 

Minimum 

0.001 

0.000 

Maximum 

0.283 

0.999 

Mean 

0.009 

0.775 

Mode 

0.008 

0.925 

Standard Deviation 

0.009 

0.199 


false color composite (RGB = Near-Infrared/Red/Green). The 
missing data in these maps correspond to 17.6 km 2 blocks 
where the MISR aerosol/surface retrievals failed because of 
cloud cover. High model-fitting RMSE has been found to 
indicate contamination by clouds [21], [22] but there is no 
indication of cloud contamination in this scene, except perhaps 
in the north-central part; it is also possible that RMSE cannot 
be used to detect cloud over snow, which also results in high 
values. The distributions of cover and height values appear 
reasonable on comparison with features in high resolution 
imagery in the Google Earth archive, i.e., high (low) canopies 
appear to have taller (shorter) trees, although this assessment 
is necessarily qualitative. 
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Fig. 12. (a) MISR fractional crown cover from inversion of the GO model versus CANAPI/QuickBird-derived fractional crown cover (b) MISR mean canopy 
heights with no rescaling versus mean LVIS RH50 canopy heights (average over the 250 x 250 m MISR cell) (c) versus RH75 (d) versus RH100. See Table II for 
inversion statistics. 



Fig. 13. (a) MISR An (nadir) camera 250 m false color composite with 
RGB = NIR, Red, Green BRFs (b) 250 m map of RMSE on GO model fitting 

(c) 250 m map of retrieved fractional crown cover, assuming crown overlap 

(d) 250 m map of retrieved mean canopy height. The dotted line indicates the 
validation area for which CANAPI and LVIS data were both available. Black 
rectangular areas indicate missing data owing to cloud. 


IV. Conclusions 

This study has demonstrated that multiangle red band surface 
BRF estimates from MISR can be used to invert a simple geo- 
metric-optical canopy reflectance model for mapping canopy 
cover and height in open forest canopies, even in areas of 
complex and often severe topography and with a very wide 
range of background conditions. The retrieved cover and height 
estimates were obtained independently from the reference data, 
not via training or empirical relationships. While a part of the 
signal observed by MISR originates from the background and is 
obtained via regression, its contribution at any given geometry 
is based on kernel weights that have a physical meaning. The 
ability of the modified simple GO model to reproduce observed 
672 nm BRF patterns was demonstrated. The use of a non-zero 
shadow component and a leaf+shadow crown signature instead 
of the non-linear Ross scattering function provided a realistic 
model, as shown by the forward-modeling results: when driven 
with first order canopy structural parameters (tree number 
density, mean crown radius, LVIS RH100 canopy heights), 
the modified model reproduced MISR BRF patterns quite 
accurately. For large area mapping there are some important 
caveats: the major factor in avoiding anomalous cover and 
height retrievals was the use of a locally-calibrated background 
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BRDF, implying that more than one set of calibration coeffi- 
cients might be required in order to enable accurate GO model 
inversion for very different landscape types. The study has also 
shown that is it difficult to obtain high precision in canopy 
height retrievals in these circumstances, with an RMSE of 6.7 m 
versus LVIS RH100. This allows only a rather coarse binning 
of canopy height to obtain high accuracy, although it is still a 
very useful capability, noting that the range of canopy heights is 
quite large (the minimum, maximum, and mean heights are 3 .4, 
62.0, and 30.1 m, respectively). Nevertheless, the study showed 
that it is possible to leverage moderate resolution multiangle 
data using a GO model to map forest structural parameters at 
regional scales; assuming data continuity, this could help to 
provide a more complete picture of changes in the forests of 
the southwestern United States as they respond to changing 
disturbance and management regimes over the coming decades. 
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